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ABSTRACT 

We measure the clustering of dark matter halos in a large set of collisionless cosmological simulations of 
the flat ACDM cosmology. Halos are identified using the spherical overdensity algorithm, which finds the 
mass around isolated peaks in the density field such that the mean density is A times the background. We 
calibrate fitting functions for the large scale bias that are adaptable to any value of A we examine. We find a 
~ 6% scatter about our best fit bias relation. Our fitting functions couple to the halo mass functions of Tinker 
et. al. (2008) such that bias of all dark matter is normalized to unity. We demonstrate that the bias of massive, 
rare halos is higher than that predicted in the modified ellipsoidal collapse model of Sheth, Mo, & Tormen 
(2001), and approaches the predictions of the spherical collapse model for the rarest halos. Halo bias results 
based on friends-of-friends halos identified with linking length 0.2 are systematically lower than for halos with 
the canonical A = 200 overdensity by ~ 10%. In contrast to our previous results on the mass function, we find 
that the universal bias function evolves very weakly with redshift, if at all. We use our numerical results, both 
for the mass function and the bias relation, to test the peak-background split model for halo bias. We find that 
the peak-background split achieves a reasonable agreement with the numerical results, but ~ 20% residuals 
remain, both at high and low masses. 

Subject headings: cosmology:theory — methods:numerical — large scale structure of the universe 



1. INTRODUCTION 

Dark matter halos are biased tracers of the underlying dark 
matter distribution. Massive halos form from high-a fluctu- 
ations in the primordial density field, inducing a correlation 
between halo mass and clus tering amp litude that is steep- 
est for cluster-sized objects (Kaiser| [l984l) . Low-mass ha- 
los are preferentially found in regions of the universe with 
below average density, thus these objects are anti-biased 
with respect to the dark matter. The clustering of galaxies 
is now underst ood through the bi as of the halos in which 
they form (e.g., Zehavi et al. 2005). Many methods that uti- 
lize galaxy clustering to constrain cosmology require precise 
knowledge of halo clustering (e.g., [v an den Bosch etalj|2003l 



Tinker et alj[2005t lAbazaiian et alj|2005t IZhe ng & WeinberL 
2007tlYoo et aTl 2009). Cosmological parameters can also be 



obtained through the abundance of high-mass halos identi- 
fied as galaxy clusters. The bias of clusters contains com- 
plementary information to their abundance. Indeed, "self- 
calibration" of cluster surveys is not possible wit hout the ad- 
ditional information present in clustering data (L ima & Hu 
l200i r2005: Majumdar & Mol u1l2004t |pgurill2009ir The pur- 
pose of this paper is to calibrate a precise, flexible halo bias 
function from numerical simulations that is accurate for dwarf 
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galaxi es through galaxy cluster masses. 

In iTinker et al .1 (1200 8) (hereafter T08), we presented a re- 
calibration of the halo mass function based on a large se- 
ries of collisionless N-body simulations. Our results uti- 
lized the spherical overdensity (SO) algorith m for identify- 
ing da rk matter halos within simulations (e.g.. lLacey & Cola 
119941) . In this approach, halos are identified as isolated den- 
sity peaks, and the mass of a halo is determined by the over- 
density A, defined here as the mean interior density relative 
to the background. Simulations of cluster formation show that 
the SO-defined halo mass should correlate tightly with clus- 
ter observable s, which are usually defined within a spherical 
apertu re (e.g.. iBialek etani200ll Ida Silva et al.1 120041 iNajall 
l2006t iKravtsov et"afll2006l) . This expectation is borne out for 
observables such as gas mass, core-excised luminosity, inte- 
grated SZ flux or its X-ray analog, Y x (e.g..lMohr et al.l|1999[ 
I Vikhlinin et al.l2006tlZhang et alJ2 008: Vik hlinin et al.120091 
Urnaud et al.ll2007ll2009HSun et al.ll2009l) . Tight correlations 
between spherical overdensity mass and observables are cru- 
cial for a robust interpretation of the observed cluster counts 
and clustering in deriving cosmological constraints. The scat- 
ter of mass-observable relations may depend on the value of 
A. In addition, particular observations may only extend out 
to a limited radius corresponding to A considerably higher 
than the often used virial value of A« 200. Thus, we seek to 
calibrate a fitting function that can be adapted to any value of 

A 

iHu & Kravtsovl (120031) and iManera et al.1 (l2009t) compared 
existing halo bias models to SO N-body results at the cluster 
mass scale. But previous studies to calibrate halo bias on nu- 
merical simulations have focused exclusivel y on the f r iends- 
of-friends (FOF) halo finding algorithm Jliiigl 1 19981 119991; 

mrienl [19991: IShefh et al.1 l200lt ISeliak & Warrenl 
" 2005L iPillepich et aiT 120081: lReedetal.1 
percolation scheme that 



Sheth & Tormen 
2004 iTinkeretalJ 
20081) . The FOF algorithm is 



makes no assumptions about halo geometry, but may spuri- 
ously group distinct halos together into the same object, con- 
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fusing the comparison between cluster observables theoreti - 
cal results dWhitefeOOlblTinker et alJ2008l:lLukic et alJ2009l) . 
Additionally, previous calibrations focus on only one value of 
the FOF linking length, / = 0.2, and thus are not applicable to 
many mass observables. Galaxy cluster studies and theoreti- 
cal halo models benefit from a self-consistently defined set of 
coupled mass and bias functions. 

The bias of halos is determined by the relative abundance of 
halos in different large-scale environments. Thus, theoretical 
models for halo bias have been derived from the mass func- 
tion using the peak- b ackground split dBardeen et ail 119861; 
Cole & Kaiser! Il98l lMo & White! 119961: IShefh & Tormenl 
19991: ISheth et ail [2001). These models produce results 
that are reasonably accurate but fail to reproduce in de- 
tail the bias of halos found in numerical simula t ions dJind 
19981: ISeliak & Warren! 12004 iTinker et alj 120051: iGao et all 
20051 iPillepich et al I I2008IT iManera et all J2009l) and 
Manera & Gaztanaga d2009l) demonstrate that using the peak 
background split to calculate the bias of massive halos from 
their mass function does not accurately match the clustering 
as measured from their spatial distribution. In addition to cal- 
ibrating the functional form of the bias, we test the peak back- 
ground split. 

In §2 we summarize our list of simulations and the numeri- 
cal techniques for calculating bias. In §3 we present our fitting 
formulae for large scale bias, comparing to previous works 
and exploring any redshift evolution. In §4 we use our results 
to test the peak-background split. In §5 we summarize our 
results. 

2. SIMULATIONS AND METHODS 

Our simulation set spans a wide range in volume and mass 
resolution in order to produce results than span nearly six 
decades in mass, from M ~ 10 10 /i"'M Q halos to massive 
clusters. The set contains 15 distinct simulations that span 
local variations of the concordance ACDM cos mology con- 
sistent with results from CMB anisotropics ( Sperg el et al.1 
l2Q03t iDunkley et al.ll2009h . Three numerical codes are rep- 
resented in our set; the adaptive refinement tech nique (ART; 
iKravtsov et al.1 1 1997 1 : iGottlober & Klypinll2008l) . the hashed 
oct-tree code (HOT; Warren et al. 2006), and th e hybrid tree- 
particle mesh code GADGET2 dSpringelll2005l) . Table 1 lists 
details of each simulation, including cosmological parame- 
ters, force resolution, volume, and mass resolution. Fur- 
ther details about the simulation set can be found in T08. 
For one simulation, L1280, there are 49 independent realiza- 
tions. The L1280 simulations were utilized in the studies on 
the halo mass func tion and bias relat i on of massive halos in 
ICrocce et all d2006l) and lManera et al.1 (|2009), as well as in the 
analysis in T08. The dark matter outputs of these simulations 
w ere kindly supplied b y R. Scoccimarro. 

ICrocce et al.1 d2006) point out that improper initial condi- 
tions can result in errors in the resulting halo mass function 
and, to a lesser extent, bias function for massive halos. These 
errors are a product of starting the simulation at too low a red- 
shift while using first-order techniques for the initial particle 
displacements and velocities. The L1280 simulations utilize 
second-order perturbation theory to ameliorate these effects. 
In T08 we performed multiple re-simulations of L1000W 
with initial conditions set using the Zel'dovich approxima- 
tion at different redshifts. We found significant differences 
between the mass functions measured in the L1000W run with 
Zel'dovich initial conditions set at z=30 and the mass function 
measured in the L1280 run. However, using Zel'dovich ini- 



tial conditions at z=60 for L1000W produces a mass function 
negligibly different from the LI 280 calculations. 

Halos are identified in each simulation using the spherical 
overdensity (SO) technique outlined in T08. In brief, the code 
identifies density peaks in the dark matter and grows spheres 
around them until the mean interior density is some set multi- 
ple, A, of the background density. Thus the mass and radius 
of a halo are related by 

M A = 3^ A p m (z)A, (1) 

where p m (z) is the mean density of the universe at redshift z. 
In our implementation of the SO algorithm, the spheres that 
contain halos are allowed to overlap; so long as the center of 
one halo is not contained within of another halo, the two 
halos are considered distinct. Owing to small overlaps in the 
exteriors of halos, a small fraction of the total mass in halos 
(~ 0.7%) is assigned to multiple halos and double counted. 
The SO method of identifying halos makes the halo mass sen- 
sitive to the force resolution of the simulation; if the density 
profile of a halo is not properly resolved, the enclosed mass 
at a given radius will be smaller. Column 10 in Table 1 lists 
the maximum value of A for which reliable results could be 
obtained for each simulation. Owing to the low spatial resolu- 
tion of the L1280 simulations, our analysis only utilizes these 
simulations for A = 200. For all simulations, only halos with 
more than 400 particles are used. This ensures that all halos 
are robustly identified and the halo profiles are well sampled. 

We define the bias of dark matter halos as the ratio of the 
halo power spectrum to the linear dark matter power spec- 
trum, 



We calculate the power spectrum of each simulation as fol- 
lows. The halos of each simulation are binned in a 200 3 den- 
sity mesh using the cloud-in-cell technique, and the power 
spectrum is computed through Fourier transformation. All 
power spectra are shot-noise subtracted. Aliasing due to the 
cloud-in-cell grid is removed through the deconvolution tech- 
nique outlined in Jing (2005). Although halo bias is scale- 
dependent in the quasi-linear and non-linear regime, here we 
focus on the large scale bias, where b is independent of k. 
We calculate b 2 as the average over the 10 largest wavelength 
modes in the simulation. For simulations with Lbox < 200 
/z -1 Mpc, non-linearity has set in at k > 10 x Itt/L^. For 
these simulations, we truncate our average to the largest 5 
modes. For the z = 2.5 outputs of two simulations, L120W and 
LI 20, the power spectra do not converge to a robust asymp- 
totic value within this fc-range, thus we exclude these outputs 
from the analysis. For each simulation, we calculate Ph(k) for 
8 jackknife subsamples of the simulation, removing density 
fluctuations from one octant of the box in each subsample. 
We use the jackknife subsamples to estimate the error on b. 

We also check these results against bias as defined by the 
halo-mass cross-power spectrum £>/,„, =Ph m /P\m- This measure 
does not require a shot noise correction, and it yields better 
statistics when the halos become vary sparse. 

We parameterize our results in terms of 'peak height' in 
the linear density field, v = 8 c /a{M), where S c is the critical 
density for collapse and a is the linear matter variance on the 
Lagrangian scale of the halo, ie, R = (3M / '4-n p m ) 1 / 3 , defined 
as 
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TABLE 1 

Properties of the Simulation Set 



Lbox h 1 Mpc 


Name 


e h 1 kpc 


N p 


m p h 1 Mq 


(Cl m ,Cl b ,cr s ,h,n) 


Code 


Zi 


Am I 


A max 


768 


H768 


25 


1024 3 


3.51 x 10 10 


(0.3,0.04,0.9,0.7,1) 


HOT 


40 





800 


384 


H384 


14 


1024 3 


4.39 x 10 9 


(0.3,0.04,0.9,0.7,1) 


HOT 


48 





3200 


271 


H271 


10 


1024 3 


1.54 x 10' 


(0.3,0.04,0.9,0.7,1) 


HOT 


51 





3200 


192 


H192 


4.9 


1024 3 


5.89 x 10 s 


(0.3,0.04,0.9,0.7,1) 


HOT 


54 





3200 


96 


H96 


1.4 


1024 3 


6.86 x 10 7 


(0.3,0.04.0.9,0.7, 1) 


HOT 


65 





3200 


1280 


L1280 


120 


640 3 


5.99 x 10" 


(0.27,0.04.0.9,0.7,1) 


GADGET2 


49 


0,0.5, 1.0 


600 


500 


L500 


15 


1024 3 x2 


8.24 X 10 9 


(0.3,0.045,0.9,0.7,1) 


GADGET2 


40 


0,0.5, 1.25, 2.5 


3200 


250 


L250 


7.6 


512 3 


9.69 x 10 9 


(0.3,0.04,0.9,0.7,1) 


ART 


49 


0, 0.5, 1.25, 2.5 


3200 


120 


L120 


1.8 


512 3 


1.07 x 10 9 


(0.3,0.04,0.9,0.7,1) 


ART 


49 


0,0.5, 1.25, 2.5 


3200 


80 


L80 


1.2 


512 3 


3.18 x 10 s 


(0.3,0.04.0.9,0.7,1) 


ART 


49 


0,0.5, 1.25, 2.5 


3200 


1000 


L1000W 


30 


1024 3 


6.98 x 10 1U 


(0.27,0.044.0.79,0.7,0.95) 


ART 


60 


0, 0.5, 1.0, 1.25 


3200 


384 


H384W 


14 


1024 3 


3.80 x 10 9 


(0.26,0.044,0.75,0.71,0.94) 


HOT 


35 





3200 


384 


H384ft„ 


14 


1024 3 


2.92 x 10 9 


(0.2,0.04,0.9,0.7,1) 


HOT 


42 





3200 


120 


L120W 


0.9 


1024 3 


1.21 x 10 8 


(0.27,0.044,0.79,0.7,0.95) 


ART 


100 


1.25,2.5 


3200 


80 


L80W 


1.2 


512 3 


2.44 X 10 8 


(0.23,0.04,0.75,0.73,0.95) 


ART 


49 


0,0.5, 1.25, 2.5 


3200 



NOTE. — The top set of 5 simulations are from the Warren et al. 12006) study. The second list of 5 simulations are of the same WMAP1 cosmology, 
but with different numerical codes. The third list of 5 simulations are of alternate cosmologies, focusing on the WMAP3 parameter set. The HOT 
code employs Plummer softening, while GADGET employs spline softening. The values of e listed for the GADGET simulations are the equivalent 
Plummer softening; when calculating the spline softening kernel, GADGET uses a value of l,4e. The force resolution of the ART code is based on the 
size of the grid cell at the highest level of refinement. A max is the highest overdensity for which halo masses can be reliably measured. 



P(k,z)W 2 (k,R)k 2 dk, 



(3) 



b{v) = 1 + ■ 



where P(k,z) is the linear power spectrum at redshift z and 
W is the Fourier transform of the top-hat window function of 
radius R. In all calculations we use 6 C = 1.686. For refer- 
ence, v of 0.75, 1, 2, and 3 (logi/ = [-0.12,0.0,0.30,0.48] H ) 
corresponds to M of 2.9 x 10 fl , 2.8 x 10 12 , 1.2 x 10 14 , and 
7.0 x 10 14 /r'M©, respectively, for the L1000W cosmology 
at z = 0. 

3. RESULTS 

3.1. Models and Measurements at A = 200 

Figure Q] shows bias as a function of v for all simulations in 
Table 1 . In this figure, halos are defined with A = 200. In the 
spherical collapse model, A « 200 defines a radius separating 
the virialized reg ion and the region of continuing infal l in an 
£l m = 1 universe (lLacev & Coldfl99H lFlelt^il l998). This 
overdensity is also close to the overdensity of halos identi- 
fied with the friends-of-friends (FOF) algorithm with the typ- 
ical linking parameter of 0.2 dDavis et al.lll98 5). Thus ana- 
lytic models are typically compared to numerical results us- 
ing either A = 200 or FOF(0.2). In Figure[T] we compare our 
A = 200 results to two current models for halo bias from the 
literature. 

First, we compare these results to predictions based on the 
spherical collapse model (SC) for the formation of dark matter 
halos. In SC, halos collapse when the linear overdensity asso- 
ciated with a peak in the densi ty field crosses a cr i tical b arrier 
5 C independent of halo mass. iPress & Schechterl (11974) used 
this model to derive an expression for the mass function of 
dark matter halos. Using the peak-backgrou nd split, which 
we wi ll de scribe in more detail in section §4, Cole & Kaiser 
(119891) and lMo & White! (1 19961) derived a bias relation of the 
form 

11 Throughout this paper, log indicates base- 10 logarithm. 



5c 



(4) 



However, the Press-Schechter mass function fails to re- 
produce the dark matter halo mass function found in simu- 
lations (see, e.g iGross et al.lll998t [Lee & Shandariiilll999l: 
I Sheth & Tormenl 119991: Uenkins et all 120011: iRobertson efall 
2009). Thus it is not surprising that the bias function in equa- 
tion dH> also does not compare well to s imulations (see, e.g., 
IImglll998LfT99llSheth & Tormenlll999 j). In Figured the SC 
model overpredicts the bias in the range 1 < v < 3, while un- 
derpredicting slightly the bias for the lowest mass halos in our 
si mulations. 

ISheth & Tormenl £[999) (hereafter ST) generalized the ex- 
pression for the Press-Schechter mass function an d calibrated 
the fre e parameters using numerical simulations. Sheth et al. 
(1200 ll) (hereafter SMT) later refined this calculation, incorpo- 
rating a "moving" barrier for the collapse criterion of halos 
in which the critical density varies with the peak height as 
motivated by the more physically realistic ellipsoidal collapse 
model. Using the peak-background split once again, SMT de- 
rived an improved expression for bias of the form 



fc(i/) = 1 + 



y/aS c 



\/a(aiy 2 ) + y/abiav 2 ) 1 



(av £ y 



2sc 



(av 2 Y + b(\-c){\-c/2)V (5) 
where a = 0.707, b = 0.5, and c = 0.6. These parameters de- 
scribe the shape of the moving barrier. In Figure [T] the SMT 
bias equation underpredicts the clustering of high-peak ha- 
los while overpredicting the asymptotic bias of low-mass ob- 
jects. The SMT function is calibrated using friends-of-friends 
(FOF) halos, thus the choice of A with which to compare is 
somewhat arbitrary, but it can be seen that the SMT function 
and our results will not agree at any overdensity: SMT bias at 
low v is too high, and is too low at high v. When increasing 
(decreasing) A, the bias at all v can only increase (decrease). 
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-0.4 -0.2 0.2 0.4 0.6 

log(i/) 

FIG. 1 . — Upper Panel: Large-scale bias as determined by the ratio (Ph/Piin) [ ^ 2 for A = 200. Results from the smaller boxes are represented by the gray circles. 
For these simulations, only measurements with less than 10% error are shown to avoid crowding. The larger-volume simulations are represented by the colored 
symbols. Each point type indicates a different simulation. The different colors, from left to right, go in order of increasing redshift from z = to z = 2.5 (see Table 
1 for the redshift outputs of each simulation). Like colors between simulations imply the same redshift. For these large-volume simulations, measurements with 
less than 25% errors are shown. Lower Panel: Fractional differences of the N-body results with the the fitting function shown in the upper panel. 



The bias formulae of ST and SM T have been shown be- 
fore to be inaccurate at low masses dSeliak & Warrenlf2 0Q4; 
iTinker et alj2005UGao et alj2005t|Piliepich et alj200oT) . Up- 
dated fitting functions have someti mes used th e functional 
form of ST dMandelbaum et all 120051) or SMT (ITinker et alj 
120051) with new parameters chosen to match numerical data, 
while others have proposed entirely new functi onal forms 
(e.g.. lSeliak & Warrenll2004l: iPillepich et al.ll2008l) . Our tests 
show that the SMT function does not yield optimal \ 2 values 
when comparing to our numerical results. We therefore intro- 
duce a similar but more flexible fitting function of the form 

b{v)=l-A^-+Bv h + Cv c . (6) 

Equation © scales as a power-law at the highest masses, fiat- 
tens out at low masses and asymptotes to b = 1 at v = 0, pro- 



vided a > 0. 

A convenient property of the SC, ST, and SMT functions is 
that they are normalized such that the mean bias of halos is 
unity. Thus, if one adopts the halo model ansatz that all mass 
is contained within halos, dark matter is not biased against it- 
self. Numerically calibrat ed bias function s in the literature do 
not obey this constraint dJindll998l 1 1999t ITinker et al.ll2005l: 
ISeliak & Warreij2004lPiilepich et al]2008l) . When fitting for 
the parameters of equation ©, we enforce this constraint by 
requiring that our bias function obey the relation 

fbMf(v)dv=l, (7) 

where f{v) is the halo mass function, once again expressed in 
terms of the scaling variable v. At each A, we use the halo 
mass functions listed in Appendix C of T08, which are nor- 
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FIG. 2. — Large-scale bias as determined by the ratio (P;, / P\ m ) 1 1 2 for four values of A. The solid line in each panel represents equation {6} with the A- 
dependent parameters listed in Table 2. The dotted curve in panel (a) is the bias formula of SMT. The dashed curve in panels (c)-(d) is the A = 200 results (i.e., 
the solid curve in panel a). 



malized such that the mean density of the universe is obtained 
when integrating over all halo masses at z = 12 . 

In T08, we found that the mass function is universal at z = 
over the range of cosmologies explored. However, the mass 
function at higher redshifts deviates systematically from the 
Z = results. In Figure [U we have included the results from 
all redshifts. Although the evolution of f(a) from z = to 
z = 1 is clear in the T08 results, the bias of these halos does not 
show significant evolution with redshift. To obtain the param- 
eters of equation (|6]), we minimize the \ 2 using the jackknife 
errors described in the previous section. The best-fit param- 
eters for the A = 200 data, listed in Table 2, yield a x 2 P er 
degree of freedom (hereafter \t) of 1.9 when incorporating 
all data from all redshifts. This high value of xt is driven 

12 The normalized mass functions in T08 are expressed in terms of \/o 
rather than v. For convenience we rewrite this function in terms of v in 
Equation {8j and give new mass function parameter values in Table 4. 



by the small error bars on the LI 280 results at z — 0. With 
49 realizations, the error bars are ~ 1% at the low -particle 
limit, thus the few percent offset between the LI 280 results 
and those of the remaining simulations yields a high xt- 
moving the LI 280 results (without refitting) yields xt = L01. 
The low spatial resolution of the L1280 simulations is a pos- 
sible source of error in the bias results. Refitting with only the 
Z = results does not change the values of the best-fit parame- 
ters or change xt- This implies that the evolution of bias with 
redshift is extremely weak between < z < 2.5, if it evolves 
at all. Simulation to simulation, the situation is not defini- 
tive. The L500 simulations shows an increased amplitude in 
the bias of — 5% between z = and z = 1 .25, but the L1000W 
simulation is consistent at all redshifts. Regardless, any evo- 
lution in the bias function at fixed v is significantly smaller 
than the evolution in the mass function. 

3.2. Large-Scale Bias as a Function of A 
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TABLE 2 
Parameters of Bias Equation 
© as a Function of A 



parameter /(A) 



A 1.0 + 0.24>'exp[-(4/v) 4 ] 

a 0.44y-0.88 

B 0.183 

b 1.5 

C 0.019+0.107y+0.19exp[-(4/y) 4 ] 

c 2.4 



NOTE. — Note: y = log 10 A. 



TABLE 3 

X?, VALUES OF THE b(l/) FITS 



A X 2 /v 



200 1.01/1.94 

300 1.33 

400 1.08 

600 1.34 

800 1.19 

1200 1.22 

1600 1.14 

2400 1.06 

3200 1.08 



NOTE. — For A = 200, the second 



value of xi includes the L1280 sim- 
ulations. 



The best-fit parameters of equation (O scale smoothly with 
A, allowing us to obtain fitting functions for these parameters 
as a function of log A. The functions that yield the parame- 
ters of equation © for 200 < A < 3200 are listed in Table 
2. Using these functions, the integral constraint in equation 
(0 is satisfied to better than 1 % for every value of A consid- 
ered. If required, higher precision can be obtained if 5 of the 
6 parameters are taken from Table 2 and the last is solved for 
numerically. 

A comparison between our numerical results and fitting 
functions for four values of A are shown in Figure|2] To avoid 
crowding and scatter, in each panel we only plot data points 
with fractional errors less than 10%. Figures [2t>-|2}i compare 
the results for A = 400, 800, and 1600 to the A = 200 fitting 
function (shown against the A = 200 data from|2^). As A 
increases, bias increases at all mass scales. At high masses 
this is expected; as A increases, a fixed set of halos will have 
lower masses but the same clustering properties, essentially 
shifting them along the i^-axis. At low masses, the amplitude 
of the bias curve also monotonically increases with A, owing 
to the substructure within high mass halos that become dis- 
tinct objects as decreases. Because these new low mass 
halos are in the vicinity of high mass objects, they have sig- 
nificant clustering. 13 

13 In principle, this makes our results sensitive to the spatial resolution of 
our simulations beyond simply resolving the halo density profiles properly. 
If subhalos are not resolved in some subset of our simulations, the change 
in bias for low u halos will be underpredicted. Our criterion for including 
simulations in our analysis is that halo density profiles are properly resolved, 
not that substructure is properly resolved. However, the fact that bias mono- 
tonically increases with A at low v is indicative that we are including these 
"revealed" subhalos. 



Table 3 shows the \ 2 V values for each value of A. The fit is 
near xt ~ 1 at all A, indicating that the fit is adequate to de- 
scribing the data even though we have combined all the simu- 
lation outputs in the fit (i.e., all cosmologies and all redshifts). 

3.3. Bias ofHigh-v Halos 

The spherical collapse model is defined by a threshold for 
collapse that is independent of halo mass. However, peaks in 
the linear density field become increasingly elliptical and pro- 
late at low v, delaying collapse. Thus, in this mass regime, the 
barrier in the ellipsoidal collapse model is significantly higher 
than the constant 8 C assumed in spherical collapse calcula- 
tions. As a result, collapsed low-mass halos reside in higher 
density environments, making them less abundant and more 
biased. At high v, the ellipsoidal collapse barrier asymp- 
totes to the spherical 5 C value and these two models should 
thus converge at high v. However, the numerically-calibrated 
barrier used in the SMT fit asymptotes to a value lower than 
the spherical collapse 5 C in order to produc e the abundance 
of hi gh-mass halos (see the discussion in Rob ertson et alj 
2009). Consequently, the clustering of high-^ halos in the 
SMT model is lower than the spherical collapse prediction. 

In Figure [3^, we compare our fitting function to the formu- 
lae of the SMT and SC models for halos with v > 1.5. We also 
show the bias results from L500 and L1000W for four differ- 
ent redshifts and from L1280 for two redshifts. We focus on 
these simulations because they are the largest in our suite 14 . 
These are the same data presented in Figure Q] but here we 
are focusing on the high-i/ regime. At v ~ 2, our simulation 
results are in good agreement with the SMT function, but at 
higher v, our results rise above the SMT function and meet 
the spherical collapse predi ction at v > 4. 

At high redshift (z ~ lOh lCohn & White! (12008) found that 
the bias of v ~ 3 halos was better described by the SC models 
rather than SMT. However, two other recent studies of halo 
bias have concluded in favor of the SMT m odel for high-peak 
halos. In contrast to lCohn & Whitel d2008l) . lR~eed et all d2008l) 
argue that the clustering of high-;/ halos at 10 < z < 30 in their 
simulations is better described by the SMT model. T hey claim 
that the bias measurements of lCohn & W hite (2008) are in er- 
ror because the bias is calculated at r = 1.5 h~ l Mpc, where the 
bias is scale-dependent. To correct for this. lReed et al.l (120081) 
use a fitting function to extrapolate the translinear correlation 
function out to linear scales. Using this techniq ue they find 
that SM T bias is a better fit than spherical collapse. iReed et aT] 
(2008) are not able to calculate error bars for their bias values, 
and the matter variance over the total volume probed in their 
simulations is ~ 12% at the redshifts for which they obtain 
their results 15 , thus s ample variance is stil l a concern. The 
numerical results of iPillepich et alJ (|2008l) are also consis- 
tent with SMT at 2 < v < 3, and deviate somewhat at higher 
masses. They use friends-of-friends (FOF) halos with a link- 
ing length of 0.2 times the mean interparticle separation, and 
they calculate halo bias by the ratio of the halo-ma tter power 
spectr um, Ph m {k), to the matter power spectrum. iGao et all 
(120051) and lAngulo et all (120081) use the halo-halo correlation 

14 We do not include the 768 h Mpc HOT simulation in this section be- 
cause the results at high v are possibly biased due to numerical issues. See 
the discussion in Appendix A of T08. We do include H768 in all fitting, but 
both the mass function and bias relation deviate from the mean results at high 
masses. 

15 The matter variance of a 1 AT 1 Mpc cube at z = 10 is 39%, but this is 
reduced by vTT due to the 1 1 realizations they have of this box. 



7 




FIG. 3. — Panel (a): The A = 200 bias function in the high-i' regime. The points with error bars represent our large- volume simulations at the redshifts listed in 
Table 1 . Only points with fractional errors less than 25% are shown. The different colors, from left to right, go in order of increasing redshift: (red, green, yellow, 
blue, cyan)=(Q.O, 0.5, 1.0, 1.25, 2.5). Like colors between simulations imply the same redshift. The dotted line is the spherical collapse prediction. The dashed 
line is the SMT function. The lower panel shows the fractional difference with respect to equation {6), 6/, = i>Nbody — ''fit ■ Panel (b): Same as (a), but now using 
bias defined by the ratio of the Phm/PlrnO 1 ). Results are shown for the L1000W simulation. Colors represent the same redshifts as in panel a. Panel (c): Bias of 
halos identified using the FOF algorithm with linking length 0.2. Bias is calculated from equation equation Results are shown fo r the L1000W simulation. 
Different colors match to different redshifts as before. The dotted curve in this Figure is the fitting function of Pillepich et al. (2008), which is calibrated on 
FOF(0.2) halos. 



TABLE 4 

Parameter of the halo mass function, equation ^ 



A 


a 
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V 


200 


0.368 


0.589 


0.864 


-0.729 


-0.243 


300 


0.363 


0.585 


0.922 


-0.789 


-0.261 


400 


0.385 


0.544 


0.987 


-0.910 


-0.261 


600 


0.389 


0.543 


1.09 


-1.05 


-0.273 


800 


0.393 


0.564 


1.20 


-1.20 


-0.278 


1200 


0.365 


0.623 


1.34 


-1.26 


-0.301 


1600 


0.379 


0.637 


1.50 


-1.45 


-0.301 


2400 


0.355 


0.673 


1.68 


-1.50 


-0.319 


3200 


0.327 


0.702 


1.81 


-1.49 


-0.336 



function to determine the bias o f FOF(0.2) halos, with results 
similar to Pill epich etail (120081) . 

In Figure^, our simulations prefer a model that is interme- 
diate between SC and SMT. But in Figure [2b we explore the 
possible systematics involved in our estimate of b{v). Here 
we use Ph m {k) to determine b(v ) from L1000W. Because shot 
noise is no longer a concern, this statistic allows us to extend 
our bias measurements to higher masses at a given redshift 
output. Although the errors at high v are large, the z = re- 
sults track our fit at all v, demonstrating that the these results 



are not due to redshift evolution (and a lack of z = data at 
v > 2). The results from other redshifts are also in agreement 
with the fit and with the z = results using Phmik). 

The last systematic to be tested is the choice of halo- 
finding algorithm. In Figure 0;, we plot the bias of halos 
in the L1000W simulation that have been identified with the 
FOF(0.2) finder. The halos were defin ed using the same algo- 
rithm and linking len gth used by both iReed et ail d2008l) and 
iPillepich et all (|2008). Although the difference with A = 200 
is small, there is a definite offset between the FOF(0.2) re- 



8 



suits and our A = 200 fit. At v ~ 3, the SMT function is a 
reasonable description of the data. At higher v the numeri- 
cal results increase faster than the v 1 scaling of SMT, but the 
errors are too large to see a sig nificant difference wit h SMT. 
The empirical fit determined by Pillepic het alJ ((2008) is also 
a good fit to our FOF(0.2) data. Their fit is consistent with 
SMT at v ~ 2 - 3 and tends higher at larger v. As discussed in 
T08 and lLukic etafl d2009b . a significant fraction of FOF ha- 
los are actually two distinct density peaks linked together by 
the FOF algorithm. This linking increases the abundance of 
massive FOF halos relative to the abundance of SO halos and 
reduces the bias. The halo b ias of the FOF(0.168) halo cat- 
alogs of Maner a"et alJ (l2009t) agree with the A = 200 results 
for the LI 280 results. 

3.4. NFW Scaling Between Values of A 

One method of obtaining halo statistics at various values of 
A is t o assume that h alos a re described by the density pro- 
file of iNavarro et alJ ( 1 19961) (hereafter NFW) and calculate 
the change in mass between the desired A and some fiducial 
valu e at which the mass fu nction or bias relation is calibrated 
(i.e.. IHu & Kravtsovll2003h . In T08 we showed that this pro- 
cedure leads to significant errors in the inferred mass function 
at M < M„, and the abundance of high mass objects is sensi- 
tive to the model for halo concentrations used. In Figure|4]we 
test this procedure on the bias function. The curves represent 
the ratio between the bias obtained using the fitting functions 
of Table 2 the bias obtained by taking the A = 200 bias func- 
tion and rescaling it to higher overdensities. We assume the 
concentration-mass relation of IZhao et alJ (l2009h for all cal- 
culations, noting that the results on the high-mass end depend 
on the model used. Models that predict a lower conc entration 
for cluster-sized halos, such as Bu llock et alJ d2001l) . yield a 
much stronger deviation from the N-body results. Scaling the 
masses from one A to another A can only result in a horizon- 
tal shift of the bias-mass relation; halos that were substructure 
at low A and revealed at high A are not taken into account. 
Low-mass substuctures in high-mass halos that are "revealed" 
as host halos when A increases will increase the mean bias of 
these objects. This is why the rescaled bias function under- 
predicts halo bias at low v. 

At the high mass end, two effects alter the agreement be- 
tween the measured bias and the rescaled bias. For the same 
object, the difference in halo mass between two values of A 
depends on the density profile of the halo. Thus, at a given 
M200, scatter in the concentration-mass relation creates a dis- 
tribution of halo masses at higher or lower A. Due to the 
steepness of the mass function at v > 1, more low bias halos 
are scattered up to higher v than high bias halos are scattered 
down. The calculation in Figure |4] assumes only the mean 
c(M) relation. Scatter accounts for half of the discrepancy. 
The rest can be accounted for by the assembly bias of halos; 
the effect that halo properties correlate with large-scale en- 
vironment (see, e.g ISheth & Torm en 2004; Gao et alj|2005t 
Idao & White! 1200ft IWechsler et al.l I20061T At v > 1, more 
concentrated halos are less clustered than on average of the 
same mass. Thus, when scaling halos of the same M200 to 
higher Ma, the value of Ma depends on c and thus depends 
on bias such that hig her values of Ma ar e less clustered. Us- 
ing the result from Wech sler et al.l (2006) that a \-a deviation 
in logc yields a ~ 13% deviation from the mean clustering 
for massive halos, scatter and assembly bias combined bring 
the rescaled bias function into agreement with the fitting func- 
tions at v > 1 . 
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FIG. 4. — The fractional difference between the bias from fitting functions 
and the bias obtained from rescaling the A = 200 fitting function to higher 
overdensities assuming NFW profiles and the concentration-mass relation of 
IZhao et alJ fiM% . 



4. TESTING THE PEAK-BACKGROUND SPLIT 

The mass function in Appendix C of T08 is written as a 
function of a. To match with our parameterization of the bias 
function in equation © and to facilitate the peak-background 
split, we rewrite this function in terms of peak height v. The 
original T08 function, g(a), is related to the new function by 
g(a) = vf(v), where 



-2A 



(8) 



Table 4 lists the values of the five parameters of equation (0 
for each value of delta. 

The mass function parameters in Table 4 are set to match 
the z = numerical results from T08. To model the redshift 
evolution of the A = 200 mass function, the parameters have 
the following redshift dependence: 



/3 = /?o(l+z)' 



0.20 



-0. (IS 



77 = 770(1+2)' 



7 = 70 (1+z)' 



0.27 



-0.01 



(9) 



(10) 



(11) 



(12) 



where (3q, etc., is the value of the parameter at z = as listed in 
Table 4. The value of a is obtained through the integral con- 
straint in equation (0. The redshift-dependent fitting function 
is accurate to ~ 5% at v > 0.6 relative to the original T08 
function. As discussed in T08, the rate of change of the mass 
function decreases as z increases, thus we recommend using 
the z = 3 in the above equations to obtain the mass function at 
z>3. 

Theoretical models for the halo mass function assume a 
one-to-one correspondence between peaks in the initial den- 
sity field and collapsed objects that form at later times. The 
peak-background split obtains the bias of halo through the 
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change in the mass function (the distribution of density peaks) 
with the large-scale density field (the background). We imple- 
ment the peak-background split under the common assump- 
tions of the excursion set formalism, such as that the smooth- 
ing mass scale (for calculating a(M) in eq uation |3l) is the 
same as the mass in the collapsed halo (see IZentnerll2007l for 
a review). 

Following ST, we define the peak height, v\, relative to the 
background, i/q, as 



[Si- 6 ] 2 



1-4 
Si 



(13) 



where on the right hand side we have only kept the leading 
order terms. We Taylor expand ^io/T^io) /vif{y{) to calculate 
the Lagrangian halo peak-background split 8\^{vi\8q). Using 
equation[8] the overabundance of halos relative to the mean in 
Lagrangian space is 



^i|<5o)i 



7 ^-(l+2ry) + 



2<p 



I + O81/1) 2 * 



— = b L (vi)5 Q . 



Hi) 



This f unction is similar to Equation 1 1 of |Sheth & T ormenl 
(1999). The Euler ian bias is related to the Lagrangian bias by 
b£ = 1 +bL- If we set Si = S c , the Eulerian bias is then 



b(v) w 1 + 



7 ^ 2 -(l+2?7) 2<\>I8 C 



l + (/3f) 20 ' 



(15) 



Figure compares peak-background split bias formula, 
equation ( fT5l ), to our N-body calibrated results using equation 
©. The peak-background split calculation does a reasonable 
job modeling the relative change in bias with A; as the nor- 
malization of the mass function is lowered by increasing A, 
the amplitude of the bias function over the mass range probed 
increases. However, at all overdensities, the peak-background 
split overestimates the bias of low-mass halos. For low over- 
densities, A < 600, equation ( TT3T > overestimates the bias of 
halos above the non-linear mass threshold. For higher over- 
densities, the N-body and analytic results appear consistent 
for high-peak halos, although the two curves must diverge 
eventually, as b ~ v 1 in the peak-background split and b ~ v 1A 
in our numerical fit. By definition, equation (TT~5T > satisfies the 
integral constraint in equation (0, as does the numerical fit; 
at log^ < -0.5, the bias from equation (|6]l is higher than the 
peak-background split calculation. 

At high masses at low overdensities, our re sults are con- 
sistent with those found in M anera et al.l (120091) . Using FOF- 
based halos, they find that employing the peak-background 
split on the mass function derived directly from their halo cat- 
alogs underpredicts the bias of high-peak halos. They find this 
result for three different values of the FOF linking parameter, 
0.2, 0.168, and 0.15. As the linking length is reduced — which 
is analogous to increasing the halo overdensity A — the dis- 
crepancy is reduced but never completely goes away. From 
T08, the linking length best associated with A = 1600 is 0.1, 
signifi cantly lower than the three values used bv lManera et alj 
(2009). Given the trends in their results and in Figure [5] we 
predict that the peak-background split will yield consistent re- 
sults at / < 0.12, but only at the high-mass end of the spec- 
trum. 

5. SUMMARY AND DISCUSSION 



We have presented a series of calibrated fitting functions 
for large-scale halo bias. The fitting functions are designed 
to yield the bias factor for any value of A within the cali- 
brated range. These functions are normalized such that, when 
used in concert with the normalized mass functions of T08 
(given by equations [8]-[12] in this paper), the overall bias of 
dark matter is unity. We find a ~ 6% scatter from simulation- 
to-simulation. Combined with the 5% error in the T08 mass 
function, this level of uncertainty has a non-negligible impact 
on the precision with which cosmological parameters can be 
constrained fro m cluster abundance studies; the Dark Energy 
Figure-of-Merit (lAlbrecht et alJl2006l) is reduced by 25-50% 
depending on th e details of the survey ( Cunh a~& Evrardl2009l : 
IWu et alj|20 09). More importantly, T08 and this study focus 
exclusively on cosmological parameters in which the vacuum 
energy density is constant with redshift. More study is re- 
quired to determine if the halo bias function is universal with 
variations in universal expansion and growth history induced 
by dark energy. For cluster studies, where the primary con- 
cern is the abundance of massive objects, a series of large- 
volume simulations similar to LIOOOW are required to ad- 
dress this uncertainty. To isolate the effects of dark energy 
in both the mass function and bias function, using the same 
initial phases with different dark energy equations of state 
would eliminate sample variance, which is a concern even for 
h~ l Gpc simulations. 

Within the precision of our data set, the numerical results 
do not show evidence for significant evolution of bias with 
redshift. Any evolution must be at the < 5% level over our 
redshift baseline. This finding contrasts with our results from 
the mass function; in T08 we demonstrated that the SO mass 
function evolves by up to ~ 50% from z = to z = 2.5. This 
evolution is more pronounced with higher overdensity. If 
the abundance of dark matter halos is connected to the bias 
of halos — as is assumed in the peak-background split — one 
would assume that b should increase at fixed v as redshift in- 
creases. To the statistical precision of our data, however, halo 
bias can be modeled by a single, redshift-independent func- 
tion. 

Although the absolute predictions of the peak-background 
split fail to reproduce our numerical results in detail, this 
method reasonably tracks the change in the bias function with 
A. Thus we can gain insight from using the peak-background 
split on the mass function at various redshifts to see how it 
changes under the peak-background ansatz. In T08, the evolu- 
tion in the mass function is mostly encompassed by a change 
in the overall normalization of vf(v) (cf., Figure 6 in T08), 
with a slightly stronger evolution for v > 1 halos. A change 
in the overall abundance of halos does not induce a change 
in their clustering. Thus, employing the peak-background 
split on the redshift-dependent mass function for A = 200 at 
z = 1-25 yields a bias function that is at nearly identical to the 
z = peak-background split function at high v, and is only 
— 5% higher at v < 0.4. 

We have paid significant attention to the bias of halos at v > 
2, which corresponds to the peak-height for galaxy clusters. 
Our A = 200 halo catalogs disfavor a bias function with an 
amplitude as low as SMT. This result is robust to any choice 
of statistic with which to ca lcula te the bias. Th e num erical 
results of lReed etall d2008l) and iPillepich et alj d2008l ) find 
good agreement with SMT at these scales, but these results 
are based on FOF(0.2) halos. The known problem of link- 
ing distinct objects in the FOF algorithm would reduce bias 
at fixed mass because two (or more) objects with intrinsically 
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FIG. 5. — Comparison of halo bias calibrated from our numerical simulations, equation j6), with results from the peak-background split, equation i\5\ . At 
A = 200, the peak-background split calculation is ~ 20% high/low and low/high v. As A increases, the residuals at v > 1 become smaller while the residuals at 
v < 1 become larger. 



lower bias are being count ed as one more ma ssive object. In 
addition, as pointed out in iLukic et all (120091) . the mean ra- 
tio between FOF halo mass and SO halo mass depends on 
concentration even for unbridged samples of halos, so this 
will also affect the relative bias between the two mass defi- 
nitions in a non-trivial manner. In our simulations, A = 200 
and FOF(0.2) do not agree. At v = 3, our FOF(0.2) results 
appear to be in agree ment with the SMT f unction as well as 
the fitting function of iPillepich et all (120081) . 

In a general sense, the peak-background split does achieve 
marked success; the first-order derivation calculated here is 
accurate to < 20% and correctly predicts the change in bias 
with A. There are several possibilities in explaining the dif- 
ferences between the theory and N-body results. For massive 
halos, a first-order expansion of t he peak-background split 
may not be sufficien t. How ever, iManera et al.1 (120091) and 
Manera & Gaztanaga (2009) demonstrate that higher-order 
terms do increase the accuracy of the calculation at high 



masses, but decreases it at lower masses. The growth of low- 
mass halos in overdensities is truncated due to the presence of 
nearby, high-mass objects dWechsler et al.l 120061 : 1 Wang et all 
I2007t iDalal et al.ll2008b lHahn et al.ll2009l) . Our implementa- 
tion of the peak-background split assumes that the local peak 
corresponding to the collapse threshold is S i = S c = 1 .686, ig- 
noring any environmental effects on the collapse of dark mat - 
ter halos. Alternatively, as discussed in Ma nera et al.l (2009), 
it is not clear that the mass that enters into the calculation of 
the peak height, S c /a(M), should be the same mass of the ob- 
ject that eventually collapses. The mass contained within the 
peak does n ot completely ma p onto the mass within the col- 
lapsed halo (IDalal et al.ll2008l) . It remains to be seen whether 
a more robust implementation of the peak-background split 
model, in which the Taylor expansion is replaced with a more 
rigorous treatment, can reconcile the differences between the- 
ory and numerical results, or if the peak-background split fails 
at a more fundamental level. More work is required to isolate 
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the failures of the model and bring our theoretical understand- 
ing of the formation of dark matter halos into agreement with 
the ever-increasing precision of numerical simulations. 
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